module nh_case2_mod
  ! Non-hydrostatic test case from Klemp, 2015
  ! Quasi Two-dimensional Mountain Ridge Case
  use const_mod

  implicit none

  private

  public nh_case2

  real(8), parameter :: Rd      = 287
  real(8), parameter :: Rv      = 461.5
  real(8), parameter :: cp      = 1004.5
  real(8), parameter :: cv      = 717.5
  real(8), parameter :: p0      = 100000
  real(8), parameter :: ztop    = 20000
  real(8), parameter :: ptop    = 3281.8_8
  real(8), parameter :: X       = 500. / 3.
  real(8), parameter :: omega   = 7.292d-5
  real(8), parameter :: lon_c   = 0.25 * pi
  real(8), parameter :: lat_c   = 0
  real(8), parameter :: h0      = 250
  real(8), parameter :: d       = 5000
  real(8), parameter :: xi      = 4000
  real(8), parameter :: peq     = 100000
  real(8), parameter :: teq     = 300
  real(8), parameter :: ueq     = 20
  real(8), parameter :: cs      = 2.5e-4
  real(8), parameter :: zh      = 10000
  real(8), parameter :: tau0    = 25
  real(8), parameter :: radius  = 6371.220d3
  real(8), parameter :: gravity = 9.80616d0

contains

  subroutine nh_case2(lon, lat, z, rho, u, v, w, pt, mr, zs, rc, g, omg, radi, zt)
    real(8), intent(in )           :: lon
    real(8), intent(in )           :: lat
    real(8), intent(in ), optional :: z
    real(8), intent(out), optional :: rho
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: w
    real(8), intent(out), optional :: pt
    real(8), intent(out), optional :: mr
    real(8), intent(out), optional :: zs
    real(8), intent(out), optional :: rc ! Rayleigh damping coefficient
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: omg
    real(8), intent(out), optional :: radi
    real(8), intent(out), optional :: zt
    
    real(8) :: dlon
    real(8) :: t, p, r, c
    
    if( ( present(rho) .or. present(u) .or. present(v) .or. present(w) .or. present(pt) .or. present(mr) .or. present(rc) ) &
    .and. .not. present(z) )then
      stop 'Need z in nh_case1 while acquiring rho, u, v, w, pt, mr, rc'
    endif
    
    if(present(zs))then
      dlon = abs( lon - lon_c )
      dlon = min( 2 * pi - dlon, dlon )
      r    = radius / x * dlon
    endif
    
    c  = 0
    
    if( present(z) )then
      t = teq
      p = peq * exp( -ueq**2 / ( 2._8 * Rd * teq ) * sin(lat)**2 - gravity * z / ( Rd * t ) )
    endif
    
    if (present(radi  )) radi  = radius / x
    if (present(rho   )) rho   = p / ( Rd * t )
    if (present(u     )) u     = ueq * cos(lat)! * sqrt( 2 * teq / t * c * z + t / teq )
    if (present(v     )) v     = 0
    if (present(w     )) w     = 0
    if (present(pt    )) pt    = t * ( p0 / p )**( Rd / Cp )
    if (present(mr    )) mr    = 0
    if (present(zs    )) zs    = h0 * exp( - r**2 / d**2 ) * cos( pi * r / xi )**2 * cos(lat)
    if (present(rc    )) rc    = merge( sin( 0.5 * pi * ( ( z - zh ) / ( ztop - zh ) )**2 ), 0._8, z>zh )
    if (present(g     )) g     = gravity
    if (present(omg   )) omg   = 0
    if (present(zt    )) zt    = ztop

  end subroutine nh_case2
end module nh_case2_mod
